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Abstract. Statistical properties of a local fluctuational fluxes measured at the 
^ ■ plasma edge are investigated in the work. It's shown that the amplitudes increments 

H , of the local fluctuational fluxes decrease by power law. For approximation of 

experimental PDFs the fractional stable distributions are used. The new algorithm 
^ . of statistical estimation of the FSD parameters based on maximum likelihood method 

• ^ I is described. By using of the algorithm the parameters of FSD are estimated by 

using experimental samples. It is shown good agreement between experimental and 
, theoretical distributions. On the basis of this results the conclusion is made about 

. applicability of the CTRW model for description of a processes underlying of the plasma 

turbulence. 
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Q^ ■ 1. Introduction 

O' 

Formation of turbulence in plasma is one of important cause which don't let to realize 
a fusion. As a investigations show, the plasma turbulence possesses have a number of 

■ properties which become apparent regardless of experiment conditions and devices. First 
^ ! of all it was established that probability density function (PDF) of amplitudes of plasma 

density fluctuation and local turbulent fluxes don't described by Gauss distribution [lH5] . 
The empirical PDFs have heavy tails with power asymptotic x"""^ and sharp peaks. 
Another property of empirical PDF is self-similarity of ones 

p{x,t) = r^p{xt-^, 1). 

This fact was established in works [2llll|6l[7] where by scaling transformation of PDFs 
was shown their identity. In works [8l|9] the Hurst parameters were calculated for the 
particle fluctuation fluxes which were measured in three stellarators TJ-IU, W7-AS, 
ATF and two tokamaks TJ-I, JET. In the works [IDIIII] have been shown that the local 
fluctuating fluxes possess self-similarity property and the increments distribution can be 
described by fractional stable distributions. In all cases, the Hurst exponent was found 
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to be if > 0.5 and lie in a fairly narrow range 0.62 ^ ii ^ 0.72 in the plasma core. 
Here H is the Hurst parameter and p(a;, 1) is distribution at the initial moment. 

The presence of the heavy-tail distributions have resulted necessity refusal of using 
classical approaches with central limit theorem and Gauss distribution in the basis 
and development new models which be able to describe the observed dependence. In 
the works [3l|T2] was proposed the model with assumption that in plasma a number of 
gaussian processes appear each of which has it own values of parameters. As a result the 
PDF of such set of gaussian processes is described by shift-scaling mixture of gaussian 
distributions. In the works the parameters of shift-scaling mixtures of distributions are 
statistically estimated with respect to experimental dates and theoretical and empirical 
distributions was compared each other. Such comparison gave a good coincidence 
theoretical and empirical distributions. 

Another method is based on proposition of cascade process dissipation of kinetic 
energy in turbulent flow. It can be explained by presence of various structures on 
different spatial scaling. This assumption leads to a multifractal conception [13l[Tl]. In 
the framework of this conception spectrum of scales can be obtained by use of structure 
function method [IHIIT^. In this work it was showed that super diffusion processes are 
developed in plasma. 

In the present work the approach with Continuous Time Random Walk (CTRW) 
model [I7j in the basis is described. This model was firstly described in the work [18] 
and further was used for description of processes in physics, chemistry, biology and 
economics. It also has found a use in the plasma physics for description of turbulent 
processes. In works [19l[20] this model lies in the basis of study of radial transport 
models of charged particles in toroidal devices. In the work [21] shown that transport 
of superthermal ions in closed magnetic trap is described by Levy walks. 

The CTRW model is based on the assumption that a randomly walking particle 
can either move or be at rest. If the PDFs of the particle jump length p{x) and the time 
during which the particle is at rest ri{t) are power-law functions 

p{x) oc ax^x'""'^, < a ^ 2, x — )• oo, 
r/(t) oc Pt^t-/^-^ < /3 ^ 1, t-^oo. 

In the framework of this assumptions the asymptotic (at t — > oo) of particle 
coordinate is described by fractional diffusion equation [17] is expressed through partial 
derivatives of fractional orders. Equations of such kind find a use for description 
of propagation of ions and electrons in the plasma of closed magnetic traps [22| [23]. 
The order of derivative with respect to time is defined by exponent /3 and order of 
derivative with respect to space coordinate is defined by exponent a. In the particular, 
for P = l,a = 2 fractional diffusion quation is reduced to the conventional diffusion 
equation 

dp{x,t) ^d'^pix,t) 
dt " dt^ ' 

where D is the diffusion coefficient. For the initial condition p{x, 0) = S{x) the solution 
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of this equation has the form 



p{x,t) = (Dt) 



/'go (x(Z^t)-^/^;2,l), 



(1) 



where gGix;2, 1) = (27r)^^/^ exp(— is the Gauss distribution. It follows from the 
([1]) that the solution is self-similar, with the Hurst exponent H = 1/2. li P ^ l,a ^ 2 
the solution of fractional diffusion equation is expressed in terms FSDs [2^ [25]. 

In the works [IDIIII] the FSDs were applied for approximation of empirical PDFs 
of local fluctuational fluxes at the edge of the plasma cord. The FSD parameters have 
been statistically estimated according to experimental time series by using methods are 
described there. Comparison of the theoretical and the empirical distributions shown 
good agreement between them. However, to estimate the distribution parameters in the 
processing of time-series had to significantly reduce the amount of the original sample. 
The reason of this is described in |10]. As a result the sample volume is about 300-500 
values, and the algorithm which was used there suppose presence a sample of an enough 
big size. As result the estimated values had enough big statistical error. In this work 
the new algorithm for statistical estimating FSD parameter is described. This algorithm 
base on maximum likelihood method and It is less demanding for the sample. By using 
this method the parameters of the FSD will estimated according to the amplitudes 
increments of the local fluctuational fluxes. After that the empirical and theoretical 
distributions will be compared together. 

The work has following structure. In the section [2] necessary theoretical information 
about FSD are presented. In the section [3] the maximum likelihood algorithm for 
statistical estimation of FSD parameter is described. In the section H] the method 
of processing of experimental time series is described and the section [5] the results of 
approximation and conclusions are presented. 

2. Fractional Stable Distributions 

Fractional stable densities (FSDs) were first introduced in p6] as limit distributions 
in the CTRW model and, independently, were used by Kolokol'tsov et al [27]. Let us 
consider the following summation scheme. Let Xi,X2, ... be independent, identically 
distributed random variables and Ti,T2, . . . be variables that are independent of both 
each other and also of the sequence Xi, X2, . . which are identically distributed over 
the positive semi-axis. We consider the compound process 



where N{t) is the counting process N{t) = T.f=iTj ^ t < Efi?^^^i> > 0. 



We assume that random variables Xj belong to the region of normal attraction of a 
strictly stable law with characteristic function g{k] a,6) {0 < a ^ 2) and Tj belong 
to the region of normal attraction of a one-sided strictly stable law with characteristic 



Nit) 




(2) 
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function g{k;(3,l) (0 < /3 ^ 1). The characteristic function of the strictly stable law 
has the form 

gik;n,e) = exp j-jfel'^exp ^-z^signA;^ | , (3) 

where < /i ^ 2, |^| ^ min(l, 2//i - 1) [28]. 

The sum ([2]) is physically interpreted as a random coordinate of a particle at time t 
that undergoes random walk in the CTRW model. Xj is then interpreted as a jump of 
the particle between two spatial points, whereas Tj is interpreted as a random waiting 
time of the particle between two successive jumps. 

In [27], it was shown that, on the above assumptions, asymptotic (at t — )■ oo) 
distribution of the sum ([2])) is described by the FSD 

oo 

g(x; a, 13,9)= [ g{xy^'^- a, e)g{y- /3, l)//"^?/, (4) 



where g[x]a,6) is the strictly stable probability density, whereas g{y;(3,l) is the one- 
sided strictly stable probability density with characteristic function ([3]). 

As already noted, the physical interpretation of the sum ([21) is the particle 
coordinate in the CTRW model. It can be shown (see [T7]) that the asymptotic 
coordinate distribution is described by the equation 

= -Di-AT^M^, t) + f^f^), (5) 

where d^/dt^ is the Riemann-Liouville fractional derivative and (— A)°/^ is the 
fractional-order Laplacian [29j. In [2l], it is shown that the solution to this equation 
has the form 

p(x,t) = {Dt^Y''\[\x\ (Z^t'^)-^/";a,/3,0) , (6) 
where q{x\ a,P, 0) is the FSD. The solution ([6]) can be represented in the form 



p{x, t) = {D't^q (\x\ (D't) ; a, /3, 

where D' = D^^^ and 

H = (3/a. (7) 

It follows herefrom that, first, the density p{x,t) possesses self-similarity and, second, 
the solution ((61) at a; — )■ oo decreases as 



3. Estimator for Parameters of FSD 



Let Yi,...,Y]\! is sample of independent identical distributed random variables. Let 
distribution of each of them coincide with FSD ([Ij) with characteristic parameters a 
and p. The task consist in estimation of a values a, (3 of the parameters a and /3 by 
sample Yi, . . . , Yjy. 
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It is well known that maximum likelihood estimation are parameters values d, /3 
under which the logarithmic likelihood function L{y; is reaching maximum: 

1 ^ 

max L(F; d, /3) = max — d, (8) 



where 



/(y,;d,/3) = logg(F,;d,/3,0) (9) 

and q{Yi] d, /3) is the FSD at the point Yi. As we can see in order to evaluate the 
maximum likelihood estimation it is necessary to be able calculate the density at the 
points Yi, . . . for given parameters values. However, in common case the FSD 
densities aren't expressed through elementary function. Therefore Monte-Carlo method 
is used for calculation of PDF of FSD. 

The task of calculation density of symmetric {9 = 0) FSD was considered in the 
work [30] and the following estimator was obtained 

where Sj{a/2) and Sj{f3) are independent identical distributed random variables either 
have one-sided strictly stable distribution with characteristic parameters a/2 and (3 
respectively. The Kanter's algorithm [31] was used for simulation of ones. 

The estimator (ITO!) allow us to estimate PDF at given point. As a result an obtained 
value don't have systematic component of an error. Moreover if we have set of points 
xi,X2, . . . , Xn, we can simultaneously estimate PDF at all given points for one realization 
of random variables Sj{a/2) and Sj{f3). It is noticeable decrease calculation volume and 
magnitude of statistical error in comparison with histogram estimation of PDF. 

Let us use the estimator (fTOl) for evaluation of likelihood function. Substitute the 
estimator in ([9]) we obtain the Monte Carlo estimator for likelihood function 



M 

i{Yi;a,/3) = -InM + ^ln 



y/AnS,{a/2) " I 45, (a/2) 



(11) 



Now the task consist in maximization of the likelihood function by parameters 

It is necessary to use direct search methods for functions of n-variables in the task 
of likelihood function maximization. In the present work the Hooke- Jeeves's method [32] 
was applied. Presence of any additional information about position of a maximum of 
function allow considerable to decrease a search time. This information will allow to set 
initial position of basic point from which the maximization algorithm will be starting. 
The nearer this point is placed toward the true value, the smaller time will necessary 
for maximization of function. It is possible to use an algorithm of statistical estimation 
of FSD parameters which was described in [33] for determination initial basic point. 



The estimators for parameters of FSD are presented in the Appendix Appendix A As 



a result a position of initial point may be determined by formulas flA.2p - flA.Sp . 
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It is necessary to take into account the fact that domain within which a maximum of 
a function is being searched is bounded. This domain is defined by domain of variation 
of parameters of FSD 

G = {(a, /3, A) : < a ^ 2, < /3 ^ 1, A > 0}. (12) 

Beyond the bounds of this domain a density of FSD isn't defined. 

It is convenient to use of helper function for searching of a maximum of likelihood 
function under conditions (fT2l) 



^^"^'^^"^ /(a;;P), P G G. ^^^^ 

Here the points P = (a, and P = (a, where a, /3, A are bounding values domain G 
for corresponding parameter. The function P) redefine likelihood function /(a:;P) 
to domain G' = {{a', /?', A') : — oo < a' < oo, — oo < < oo, — oo < A' < oo} and it is 
penalty for exit of algorithm out of the domain G. This function is decreasing the value 
of the function l{x] P) if point P ^ G. At the same time such decreasing the stronger, 
the farther the point P exceed the bounds of the domain G. If the point P G G then 
l{x; P) is not changing. Thus, the function (fT3|l decrease a value of likelihood function if 
the algorithm exceed the bounds of the G domain. This not allow to algorithm to move 
further in this direction. As a result this lead to returning the algorithm back into the 
G domain. 



4. Approximation of turbulent particle fluxes 

Now we apply the FSD to describing probability densities of amplitudes of turbulent 
particle fluxes. 

Turbulent particle fluxes were measured in the edge plasma of the L-2M stellarator 
[M] with the use of a set of three Langmuir probes. It is well known that the flux F^- is 
defined as, 

Tj = 6nj6vj, j = 1 . . . N, 

where 6nj are plasma density fluctuations, 6v^ are fluctuations of the radial velocity, 
and is the number of points in the time sample being processed. The fluctuations of 
the radial velocity are computed by the formula Sv'j = c{6ip^ — 6(p'j)/rBA9, where 6ip^ 
and 6(p'^ are the floating-potential fluctuations measured by the first and second probes, 
r is the mean radius of the magnetic flux surface, c is the speed of light, and B is the 
magnetic field. 

Hence, the problem reduces to finding estimates a, /3, A of parameters of the 
FSD from measured time sequences A by using the algorithm (JHl). However, it 
is incorrect to apply this algorithm to time sequences Tj. First, we deal with real 
physical variables that vary in time at a random but finite rate. The finite rate 
implies that two successive variables in the measure time sample are not statistically 
independent. Second, an analog-to-digital converter (ADC) introduces additive noise 
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into the measurement results. This noise disturbs homogeneity of the sample and, 
consequently, we cannot use the algorithm for estimation of the FSD parameters. We 
have to sample statistically independent variables and to eliminate the noise from the 
time sequences. 

We assume that the observed plasma turbulence is due to the presence of structures 
that have an enhanced density of charged particles in comparison with the surrounding 
plasma. We assume that the density of charged particles in a given structure does not 
depend on the density in the neighboring structures. As a result at the instant at which 
the structure passes near the probe, its signal is at a maximum, and it is at a minimum 
when the structure leaves this region. Choosing only maximum and minimum values of 
the probe current and passing to their increments, we obtain a sample of independent, 
identically distributed random variables AF'^ = F'^^ — F'^^_j^, 

It has been found that the ADC noise does not exceed a certain maximal Cmax 
and minimal Cmin levels, which can be determined from a signal measured in the 
absence of a plasma (before the ECR heating pulse is switched on). The procedure 
is as follows: before the heating pulse, we determine values Cmax and Cmin and compute 
AC = I Cmax — Cmin I- From the sample AF'^ we exclude increments that satisfy the 
condition | AF'| < AC, and proceed to the new sample AF^, m = 1, . . . , A^": 

Af; = Af'3',_2 - 1/2 (Af'3',_, + Af:;,) , (14) 

where / = 1, . . . , [A^"/3]. Here, [A\ is the integer part of A (for more details see |10j). 

From the sample AF and formulas (|8]) we can find estimates a, /3, A of the 
parameters a, /3, A. 

5. Results and Conclusions 

The method is described above for approximation of local fluctuational fluxes was used. 
The fluxes were obtained in the set of experiments [31] on the device L-2M. The three 
parameters Q!,/3 and A were estimated. The results of approximation are presented in 
the Table [H Also on the Figure [1] show the empirical probability densities of increments 
of amplitudes of turbulent particle fluxes AF*. It should be noted that these figures 
show the increments of amplitudes computed by formula dHD, instead of the increments 
of amplitudes of physical variables under study. It can be seen from these curves that 
the empirical probability density distributions have tails decreasing as which is 

consistent with the results of other authors [H [2|f¥P6U8|JT2] . Comparing these distributions 
with FSDs shows that, within the statistical errors, the FSDs give good fit to the 
empirical distributions. 

The x^-criterion was used to test the hypothesis about coincidence of theoretical 
and empirical PDFs. In all cases the following value of the significance level p = 0.01 
and the number of freedom degrees m = 10 were used. Theoretical value is Xt = 23.21 in 
all cases. Calculated values of xl magnitude in each cases are presented in the Table [H 
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Table 1. The estimation results for fluctuation fluxes which were carried out in the 

work are presented in the table. 

No. d /? A xl No. d /3 A xl 

50658 f.54 0.29 0.70 f6.82 53080 1.28 0.42 0.20 9.31 

50659 1.58 0.20 0.70 22.07 53125 1.34 0.30 0.2 17.47 
50661 1.26 0.50 0.40 10.61 53129 1.18 0.62 0.12 11.99 
50664 1.31 0.19 0.50 15.00 53134 1.38 0.03 0.20 9.80 




Figure 1. The PDFs of increments of fluctuational fluxes are shown. On the top 
figure the discharge No. 50659 is shown and on the bottom figure discharge No. 

53080 is shown. • are empirical PDF and are FSDs with estimated values of the 

parameters (see Table [T]). 

It can be seen that in all cases Xe < X? • This means that the results don't contradict 
the hypothesis of the possibility of approximation of local fluctuation fluxes by FSDs. 

The FSD, it will be remembered, is a solution to the generalized diffusion equation 
(jS]), hence the coincidence of FSD with the empirical density means that the increments 
of fluctuation amplitudes obey the generalized diffusion equation. This allows the 
following conclusions. As is well known, equation ([5]) describes the asymptotic 
distribution of the particle coordinate in the CTRW model. This means that stochastic 
process underlying the turbulent particle fluxes can formally be described in the CTRW 
model. It should be borne in mind that the walk occurs in the phase space (the time 
is a characteristic of interest). We use the space (t, F) for Figure [H Accordingly, the 
coordinate x in equation ^ will take physical variables F. 
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Appendix A. Estimation of the parameters by moment method 

Let Zi, ^2, . . . , Zni n ^ 4 be independent, identicaUy distributed random variables with 
density dl]). The problem is to determine estimates A of unknown parameters 

a, 13,6, X. This problem was solved in [33], where a factional stable stochastic variable 
was represented in the form 

where Y{a, 6) and S'(/3, 1) are strictly stable (0 < a ^ 2) and one-sided strictly stable 
(0 < /3 ^ 1) stochastic variables with characteristic function ([3]). 

Here, we only present the final result. The formulas for estimates a, /3, 9, X of 
parameters a, /3,9, X has the form 

2 " 

^ =1--5^I(X,<0), (A.2) 

27r , 
a = (A.3) 



12K + (2Zn + 39^-1 

/3 = (A.4) 
A=exp{f/„-C(Z„-l)}, (A.5) 

/ \ 1/3 

where Zn = f 1 + 2^^) , Un, Vn, are sample centered logarithmic moments 



^ n 1 " 

f/„, =-5^1n|X,|,K = -5^(ln|X,|-f//^ 

n 

-5^(ln|X,|-t/, 



n) 1 



n 



3 

n) 1 



1(A) is the indicator of event A, (D = 0.5772156649015325 ... is the Eulerian constant, 
and C(3) is the Riemann function at point 3. 
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